home *** CD-ROM | disk | FTP | other *** search
/ Libris Britannia 4 / science library(b).zip / science library(b) / DDJMAG / DDJ9203.ZIP / 80X87.ZIP / QUAD.ASM < prev    next >
Assembly Source File  |  1992-01-24  |  3KB  |  81 lines
  1.  
  2. ; quad.asm: integer power function callable from Borland C++.
  3. ; Copyright (C) 1991 by Nicholas Wilt.  All rights reserved.
  4.  
  5. .MODEL  LARGE,C
  6.  
  7. .CODE
  8.  
  9. ; int solve_quadratic(double a, double b, double c, double *x1, double *x2);
  10. ; solve_quadratic takes the coefficients of a quadratic polynomial
  11. ; and finds the roots of that polynomial.  If there are two real
  12. ; roots, it writes them back to x1 and x2 and returns 1.  If the
  13. ; discriminant b^2-4*a*c is less than 0, the roots are not real
  14. ; and solve_quadratic returns 0.
  15. ; The quadratic formula is as follows:
  16. ;   (-b +/- sqrt(b*b-4*a*c)) / 2*a
  17.  
  18.     PUBLIC  solve_quadratic
  19.  
  20. solve_quadratic PROC    A:QWORD,B:QWORD,C:QWORD,X1:DWORD,X2:DWORD
  21.                 ; Comments show stack contents
  22.                 ; separated by | signs
  23.                 ; Stack top is at left
  24.     fld A       ; a     Here, ST(0) contains a.
  25.     fld B       ; b | a     Now ST(1) has a.  Etc.
  26.  
  27.     ; Negate b -- squaring it is sign-independent, and we need b
  28.     ; negated in all the other instances in the formula.
  29.     fchs
  30.  
  31.     fld C       ; c | b | a
  32.     fld st(1)       ; b | c | b | a
  33.     fmul    st,st(0)    ; b*b | c | b | a
  34.  
  35.     ; Just plain fxch has an implicit operand of ST(1).
  36.  
  37.     fxch            ; c | b*b | b | a
  38.     fadd    st,st(0)    ; 2*c | b*b | b | a
  39.     fadd    st,st(0)    ; 4*c | b*b | b | a
  40.     fmul    st,st(3)    ; 4*a*c | b*b | b | a
  41.     fsubp   st(1), st   ; b*b-4*a*c | b | a
  42.     ftst            ; Compare against 0
  43.  
  44.     ; We've computed b*b-4*a*c.  If negative, we return 0.
  45.     ; To do the comparison, we have to store the 80x87 status
  46.     ; word and use sahf to store it into the flags.  Once it's
  47.     ; in the flags, we can jump on above or equal (jae) to jump if the
  48.     ; number tested is greater than 0 after the FTST instruction above.
  49.     sub sp,2        ; Quick, allocate a local
  50.     mov bx,sp       ; Point BX at it
  51.     fstsw   ss:[bx]     ; Store the 80x87's status word there
  52.     mov ax,ss:[bx]  ; AX <- status word
  53.     add sp,2        ; Deallocate the local
  54.     sahf            ; Get SW into flags
  55.     jae ComputeResults  ; Jump if positive
  56.     fstp    st      ; This instruction clears the stack
  57.     fstp    st      ; we've got three values on the stack
  58.     fstp    st      ; to clear
  59.     xor ax,ax       ; Return 0 - no roots found.
  60.     jmp short LeaveQuadratic
  61. ComputeResults:
  62.     fsqrt           ; Find square root of discriminant
  63.     fxch    st(2)       ; a | b | sqrt(b*b-4*a*c)
  64.     fadd    st,st(0)    ; 2*a | b | sqrt(b*b-4*a*c)
  65.     fld st(1)       ; b | 2*a | b | sqrt(b*b-4*a*c)
  66.     fadd    st,st(3)    ; b+sqrt(b*b-4*a*c) | 2*a | b | sq..
  67.     fdiv    st,st(1)    ; x1 | 2*a | b | sqrt(b*b-4*a*c)
  68.     les bx,X1       ;
  69.     fstp    qword ptr es:[bx] ; 2*a | b | sqrt(b*b-4*a*c)
  70.     fxch            ; b | 2*a | sqrt(b*b-4*a*c)
  71.     fsubrp  st(2),st    ; 2*a | -b - sqrt(b*b-4*a*c)
  72.     fdiv            ; x2
  73.     les bx,X2       ; Store x2
  74.     fstp    qword ptr es:[bx]
  75.     mov ax,1        ; Return 1
  76. LeaveQuadratic:
  77.     ret         ; Return
  78. solve_quadratic ENDP
  79.  
  80.     END
  81.